反応拡散系モデルのアニメーション¶

なんか生物系のモデルでおもしろそうな絵が描けるモデルがあったので実装した。

反応拡散系のモデルで、Gray-Scottモデルとかいうらしい(Gray and Scott 1983, 1984)。

\begin{align} \frac{\partial A}{\partial t}&=D_{A} \Delta A -AB^2 + F(1-A) \label{A}\tag{1} \\ \frac{\partial B}{\partial t}&=D_{B} \Delta B +AB^2 -(F+k) B \label{B}\tag{2} \end{align}

この式は、二つの物質が反応しながら拡散していくモデルらしい。その反応は、簡単に

  1. 物質Aが供給
  2. 2つのBが1つのAを使って3つのBになる：$A+2B \rightarrow 3B$
  3. 物質Bが取り除かれる

という感じらしい。右辺第1項は拡散、第2項は$A+2B \rightarrow 3B$の反応を表し、(1)式の右辺第3項によって供給され、(2)式の右辺第3項によって取り除かれる。拡散係数$D_A,D_B$と供給・減量を表す係数$F,k$によってパラメタライズされる。

このモデルを、素朴に差分法で実装した。周期境界条件を課しており、ラプラシアンは周り9点で計算した。

あんまり重くなるのも嫌なので、グリッド数は100$\times$100にして、5000回時間発展させた。

かなり愚直に書いているので、もっと高速化する余地がありそう。

In [1]:
using Plots
In [2]:
function main(; DA=1.0, DB=0.5, f=0.0545, k=0.062)
#     DA=1.0
#     DB=0.5
#     f=0.0545
#     k=0.062
    dt=1.0

    nx=100
    ny=100
    N=5000
    freq=100

    A=ones(nx+2,ny+2)
    B=zeros(nx+2,ny+2)
    r=20
    A[nx÷2-r:nx÷2+r, ny÷2-r:ny÷2+r].=0.5
    B[nx÷2-r:nx÷2+r, ny÷2-r:ny÷2+r].=0.25
    A+=0.05*rand(nx+2,ny+2)
    B+=0.05*rand(nx+2,ny+2)
    A[1   ,:]=A[nx+1,:]
    A[nx+2,:]=A[2   ,:]
    A[:,1   ]=A[:,ny+1]
    A[:,ny+2]=A[:   ,2]
    B[1   ,:]=B[nx+1,:]
    B[nx+2,:]=B[2   ,:]
    B[:,1   ]=B[:,ny+1]
    B[:,ny+2]=B[:   ,2]

    A_sub=zeros(nx+2,ny+2)
    B_sub=zeros(nx+2,ny+2)
    
    As = Array{Array{Float64,2}}(undef, N÷freq+1)
    Bs = Array{Array{Float64,2}}(undef, N÷freq+1)
    
    num=1
    As[num]=copy(A)
    Bs[num]=copy(B)

    for step=1:N
        ABB=A.*B.*B
        A_sub[2:nx+1,2:ny+1]=(DA*(
                -A[2:nx+1,2:ny+1]+
                0.2*(A[2:nx+1,3:ny+2]+A[2:nx+1,1:ny]+A[3:nx+2,2:ny+1]+A[1:nx,2:ny+1])+
                0.05*(A[3:nx+2,3:ny+2]+A[3:nx+2,1:ny]+A[1:nx,3:ny+2]+A[1:nx,1:ny])
                ) - ABB[2:nx+1,2:ny+1] + f*(1.0.-A[2:nx+1,2:ny+1]))*dt
        B_sub[2:nx+1,2:ny+1]=(DB*(
                -B[2:nx+1,2:ny+1]+
                0.2*(B[2:nx+1,3:ny+2]+B[2:nx+1,1:ny]+B[3:nx+2,2:ny+1]+B[1:nx,2:ny+1])+
                0.05*(B[3:nx+2,3:ny+2]+B[3:nx+2,1:ny]+B[1:nx,3:ny+2]+B[1:nx,1:ny])
                ) + ABB[2:nx+1,2:ny+1] - (k+f)*B[2:nx+1,2:ny+1])*dt
        A_sub[1   ,:]=A_sub[nx+1,:]
        A_sub[nx+2,:]=A_sub[2   ,:]
        A_sub[:,1   ]=A_sub[:,ny+1]
        A_sub[:,ny+2]=A_sub[:   ,2]
        B_sub[1   ,:]=B_sub[nx+1,:]
        B_sub[nx+2,:]=B_sub[2   ,:]
        B_sub[:,1   ]=B_sub[:,ny+1]
        B_sub[:,ny+2]=B_sub[:   ,2]

        A=A+A_sub
        B=B+B_sub
        
        if step%freq==0
            num+=1
            As[num] = copy(A)
            Bs[num] = copy(B)
            print(step, " ")
#             println(A[2,2])
        end
    end
#     ha=heatmap(A)
#     hb=heatmap(B)
#     plot(ha, hb, layout=(1,2), size=(1000,400))
    return As,Bs
end
Out[2]:
main (generic function with 1 method)
In [3]:
@time As,Bs=main();

100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000  10.555633 seconds (557.62 k allocations: 17.982 GiB, 16.07% gc time, 0.39% compilation time)
In [4]:
N=length(As)
anim = @animate for t=1:N
    ha=heatmap(As[t],clim=(0.5,1),c=:matter,title="A")
    hb=heatmap(Bs[t],clim=(0,0.5),c=:matter,title="B")
    plot(ha, hb, layout=(1,2), size=(1000,400))
end
filename="reaction.gif"
gif(anim, filename, fps=10)

[ Info: Saved animation to C:\****\****\****\****\reaction.gif
Out[4]:

なんかウネウネしてて面白い。パラメータを変えればほかの挙動も見える。

周期境界条件ゆえにフーリエ変換して波数空間で時間発展させてもいいかもしれないので、そっちも余裕があれば書きたい。

参考¶

Reaction-Diffusion Tutorial